22/08/24
Quick revisit of old ePS results (2019) to compare with ELI measurements.
Working from recent N2O notebooks, but may still need some plotter updates. (See https://phockett.github.io/ePSdata/N2O-preliminary/N2O_orbs6-9_LF-MF_preliminary_220224-tidy.html.)
For pol state/rotations, see https://epsproc.readthedocs.io/en/3d-afpad-dev/special_topics/frame_definitions_docs_130323.html#Pol-with-the-Epol-class and https://phockett.github.io/ePSdata/N2O-preliminary/N2O_orbs6-9_MF-polStates_180324.html.
20/03/24
Added basic MF-PAD results (linear polarization, $(x,y,z)$ alignments).
22/02/24 PH
1st processing for N2O results.
For methods: https://epsproc.readthedocs.io/en/dev/demos/ePSproc_class_demo_161020.html
%matplotlib inline
# Quick hack to override default HTML template
# NOT required in some JLab versions.
# https://stackoverflow.com/a/63777508
# https://stackoverflow.com/questions/21971449/how-do-i-increase-the-cell-width-of-the-jupyter-ipython-notebook-in-my-browser
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))
/tmp/ipykernel_80233/843188395.py:7: DeprecationWarning: Importing display from IPython.core.display is deprecated since IPython 7.14, please import from IPython display from IPython.core.display import display, HTML
import sys
import os
from pathlib import Path
import numpy as np
# import epsproc as ep
import xarray as xr
import matplotlib.pyplot as plt
from datetime import datetime as dt
timeString = dt.now()
import epsproc as ep
# Plotters
from epsproc.plot import hvPlotters
# Multijob class dev code
from epsproc.classes.multiJob import ePSmultiJob
hvPlotters.setPlotters(width = 700, snsStyle='whitegrid')
OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.
* sparse not found, sparse matrix forms not available. * natsort not found, some sorting functions not available.
* Setting plotter defaults with epsproc.basicPlotters.setPlotters(). Run directly to modify, or change options in local env.
* Set Holoviews with bokeh. * pyevtk not found, VTK export not available.
* Set Holoviews with bokeh.
# For class, above settings don't take, not sure why, something to do with namespaces/calling sequence?
# Overriding snsStyle does work however... although NOT CONSISTENTLY????
# AH, looks like ordering matters - set_style LAST (.set seems to override)
import seaborn as sns
sns.set(rc={'figure.figsize':(10,6)}) # Set figure size in inches
sns.set_context("paper")
sns.set_style("whitegrid") # Set plot style
sns.set_palette("Paired") # Set colour mapping
# Try direct fig type setting for PDF output figs
from IPython.display import set_matplotlib_formats
# set_matplotlib_formats('png', 'pdf')
set_matplotlib_formats('svg', 'pdf')
/tmp/ipykernel_80233/1412183033.py:14: DeprecationWarning: `set_matplotlib_formats` is deprecated since IPython 7.23, directly use `matplotlib_inline.backend_inline.set_matplotlib_formats()`
set_matplotlib_formats('svg', 'pdf')
# xr.set_options(display_style='html')
# import warnings
# # warnings.filterwarnings('once') # Skip repeated numpy deprecation warnings in current build (xr15 env)
# warnings.filterwarnings('ignore') # Skip repeated numpy deprecation warnings in current build (xr15 env)
# # Scan for subdirs, based on existing routine in getFiles()
# fileBase = Path('/home/paul/ePS/OCS/OCS_survey') # Data dir on Stimpy
# fileBase = Path('/home/paul/fock-mount/globalhome/eps/N2O/N2O_valence') # Data dir on Jake
fileBase = Path('/home/jovyan/jake-home/femtobackSSHFS/FemtoShare/Results/ePolyScat_stuff/co2') # Data dir on Jake via Docker QM3 container
# TODO: fix orb label here, currently relies on (different) fixed format
data = ePSmultiJob(fileBase, verbose = 0, ext='.inp.out', jobStructure='subDirs') # TODO: fix subDirs usage!
# CO2 case: fix by setting .inp.out, otherwise also reads base dir (older results)
data.scanFiles(outputKeyType='orb')
data.jobsSummary()
Found 4 directories, with 4 files.
*** Job orb10 details
Key: orb10
Dir /home/jovyan/jake-home/femtobackSSHFS/FemtoShare/Results/ePolyScat_stuff/co2/co2_1pg_0.1-50.1eV, 1 file(s).
{ 'batch': 'ePS co2, batch co2_1pg_0.1-50.1eV, orbital E1G',
'event': ' CO2 1pig-1',
'orbE': -14.871022583432442,
'orbLabel': 'orb10',
'symmetryLabel': '# CO2 1pig-1'}
*** Job orb8 details
Key: orb8
Dir /home/jovyan/jake-home/femtobackSSHFS/FemtoShare/Results/ePolyScat_stuff/co2/co2_1pu_0.1-50.1eV, 1 file(s).
{ 'batch': 'ePS co2, batch co2_1pu_0.1-50.1eV, orbital E1U',
'event': ' CO2 1piu-1',
'orbE': -19.79084121670707,
'orbLabel': 'orb8',
'symmetryLabel': '# CO2 1piu-1'}
*** Job orb7 details
Key: orb7
Dir /home/jovyan/jake-home/femtobackSSHFS/FemtoShare/Results/ePolyScat_stuff/co2/co2_3sigu_0.1-50.1eV, 1 file(s).
{ 'batch': 'ePS co2_2016, batch co2_3sigu_0.1-50.1eV, orbital A2U',
'event': ' CO2 3sigu-1',
'orbE': -20.35955918924822,
'orbLabel': 'orb7',
'symmetryLabel': '# CO2 3sigu-1'}
*** Job orb6 details
Key: orb6
Dir /home/jovyan/jake-home/femtobackSSHFS/FemtoShare/Results/ePolyScat_stuff/co2/co2_4sigg_0.1-50.1eV, 1 file(s).
{ 'batch': 'ePS co2, batch co2_4sigg_0.1-50.1eV, orbital A1G',
'event': ' CO2 3sigu-1',
'orbE': -21.709243947049224,
'orbLabel': 'orb6',
'symmetryLabel': '# CO2 3sigu-1'}
*** Job stacked details
Key: stacked
No 'job' info set for self.data[stacked].
# # Fix labels for plots
# for key in data.data.keys():
# data.data[key]['jobNotes']['orbLabel'] = key
Note orbital numbering in table below:
Orb is Gamess file output orbital numbering, ordered by energy but not grouped by degeneracy.OrbGrp is grouped numbering by degeneracy, also used by ePolyScat, and will be used for labels etc. below.data.molSummary()
*** Molecular structure
*** Molecular orbital list (from ePS output file) EH = Energy (Hartrees), E = Energy (eV), NOrbGrp, OrbGrp, GrpDegen = degeneracies and corresponding orbital numbering by group in ePS, NormInt = single centre expansion convergence (should be ~1.0).
| props | Sym | EH | Occ | E | NOrbGrp | OrbGrp | GrpDegen | NormInt |
|---|---|---|---|---|---|---|---|---|
| orb | ||||||||
| 1 | SU | -20.6434 | 2.0 | -561.735531 | 1.0 | 1.0 | 1.0 | 0.974763 |
| 2 | SG | -20.6433 | 2.0 | -561.732810 | 1.0 | 2.0 | 1.0 | 0.977603 |
| 3 | SG | -11.4445 | 2.0 | -311.420710 | 1.0 | 3.0 | 1.0 | 1.000000 |
| 4 | SG | -1.5454 | 2.0 | -42.052476 | 1.0 | 4.0 | 1.0 | 0.999082 |
| 5 | SU | -1.4923 | 2.0 | -40.607552 | 1.0 | 5.0 | 1.0 | 0.998756 |
| 6 | SG | -0.7978 | 2.0 | -21.709244 | 1.0 | 6.0 | 1.0 | 0.999632 |
| 7 | SU | -0.7482 | 2.0 | -20.359559 | 1.0 | 7.0 | 1.0 | 0.999765 |
| 8 | PU | -0.7273 | 2.0 | -19.790841 | 1.0 | 8.0 | 2.0 | 0.999973 |
| 9 | PU | -0.7273 | 2.0 | -19.790841 | 2.0 | 8.0 | 2.0 | 0.999973 |
| 10 | PG | -0.5465 | 2.0 | -14.871023 | 1.0 | 9.0 | 2.0 | 0.999965 |
| 11 | PG | -0.5465 | 2.0 | -14.871023 | 2.0 | 9.0 | 2.0 | 0.999965 |
*** Warning: some orbital convergences outside single-center expansion convergence tolerance (0.01): [[1. 0.97476329] [2. 0.97760302]]
TODO:
# With defaults
# eDiag = ep.setPolGeoms(defaultMap='exyDiag') # Default mapping
# data.MFBLM(RX = eDiag)
# Test alternative z mappings - custom case
tPoints = 10
tRot = np.linspace(0,np.pi/2,tPoints)
pRot = np.zeros(tPoints)
cRot = np.zeros(tPoints)
# labels = tRot.round(2).astype(str)
labels = tRot.round(2) #.astype(str)
eulerAngs = np.array([labels, pRot, tRot, cRot]).T
eRotRX = ep.setPolGeoms(eulerAngs)
# pRot = [0, np.pi/4]
# tRot = [0, np.pi/2]
# cRot = [0, 0]
# labels = ['x', 'y', 'z']
# Single orb
# orbKeys = 'orb10'
# data.MFBLM(keys=orbKeys, RX = eRotRX)
# All orbs
data.MFBLM(RX = eRotRX)
# Plot all cases with widgets
data.BLMplot(dataType='MFBLM', thres = 1e-2, backend='hv', width=700, ylim=(-1.5, 2.5), xlim=(0,15))
Applying MFBLM plot defaults, pass overrideMF=True to skip. BLMplot set data and plots to self.plots['BLMplot']
True
# Plot single Eke vs. pol rotation ("label" coord)
data.BLMplot(dataType='MFBLM', thres = 1e-2, backend='hv', width=700, ylim=(-1.5, 2.5), xDim='Labels', # Erange=[0,2],
selDims={'Eke':[10.1]}) #, Etype=None) #,keys='orb10')
Applying MFBLM plot defaults, pass overrideMF=True to skip. BLMplot set data and plots to self.plots['BLMplot']
True
Update 30/08/24: recalc for both degen components.
Plots show separate components, then sum over components.
orbKey = 'orb10'
# Recalc for it=1,2 components
data.MFBLM(keys=orbKey, selDims={'it':[1,2],'Type':'L'},RX = eRotRX)
data.padPlot(keys=orbKey, Erange = [10, 15, 5], selDims={'it':1, 'Labels':tRot.round(2)[[0,2,5,7,9]]}, dataType='MFBLM', backend='pl', returnFlag=True,
height=400, width=1200, rc=[1,5]) #, facetDims=['Eke','Labels']) #['Eke','Labels'])
data.padPlot(keys=orbKey, Erange = [10, 15, 5], selDims={'it':2, 'Labels':tRot.round(2)[[0,2,5,7,9]]}, dataType='MFBLM', backend='pl', returnFlag=True,
height=400, width=1200, rc=[1,5]) #, facetDims=['Eke','Labels']) #['Eke','Labels'])
Using default sph betas. Summing over dims: set() Plotting from self.data[orb10][MFBLM], facetDims=['Eke', 'Labels'], pType=a with backend=pl.
Set plot to self.data['orb10']['plots']['MFBLM']['polar'] Using default sph betas. Summing over dims: set() Plotting from self.data[orb10][MFBLM], facetDims=['Eke', 'Labels'], pType=a with backend=pl. *** WARNING: plot dataset has min value < 0, min = (-0.00032201299064423396-1.9702018625549438e-17j). This may be unphysical and/or result in plotting issues.
Set plot to self.data['orb10']['plots']['MFBLM']['polar']
# With sum over degen components
data.padPlot(keys=orbKey, Erange = [10, 15, 5], selDims={'Labels':tRot.round(2)[[0,2,5,7,9]]}, dataType='MFBLM', backend='pl', returnFlag=True,
height=400, width=1200, rc=[1,5]) #, facetDims=['Eke','Labels']) #['Eke','Labels'])
Using default sph betas.
Summing over dims: {'it'}
Plotting from self.data[orb10][MFBLM], facetDims=['Eke', 'Labels'], pType=a with backend=pl.
Set plot to self.data['orb10']['plots']['MFBLM']['polar']
Update 30/08/24: recalc for both degen components.
Plots show separate components, then sum over components.
orbKey = 'orb8'
# Recalc for it=1,2 components
data.MFBLM(keys=orbKey, selDims={'it':[1,2],'Type':'L'},RX = eRotRX)
data.padPlot(keys=orbKey, Erange = [10, 15, 5], selDims={'it':1, 'Labels':tRot.round(2)[[0,2,5,7,9]]}, dataType='MFBLM', backend='pl', returnFlag=True,
height=400, width=1200, rc=[1,5]) #, facetDims=['Eke','Labels']) #['Eke','Labels'])
data.padPlot(keys=orbKey, Erange = [10, 15, 5], selDims={'it':2, 'Labels':tRot.round(2)[[0,2,5,7,9]]}, dataType='MFBLM', backend='pl', returnFlag=True,
height=400, width=1200, rc=[1,5]) #, facetDims=['Eke','Labels']) #['Eke','Labels'])
Using default sph betas. Summing over dims: set() Plotting from self.data[orb8][MFBLM], facetDims=['Eke', 'Labels'], pType=a with backend=pl. *** WARNING: plot dataset has min value < 0, min = (-0.0027680553382225343+5.606882925141018e-17j). This may be unphysical and/or result in plotting issues.
Set plot to self.data['orb8']['plots']['MFBLM']['polar'] Using default sph betas. Summing over dims: set() Plotting from self.data[orb8][MFBLM], facetDims=['Eke', 'Labels'], pType=a with backend=pl. *** WARNING: plot dataset has min value < 0, min = (-0.0015690310397590967-1.5395670849294163e-17j). This may be unphysical and/or result in plotting issues.
Set plot to self.data['orb8']['plots']['MFBLM']['polar']
# With sum over degen components
data.padPlot(keys=orbKey, Erange = [10, 15, 5], selDims={'Labels':tRot.round(2)[[0,2,5,7,9]]}, dataType='MFBLM', backend='pl', returnFlag=True,
height=400, width=1200, rc=[1,5]) #, facetDims=['Eke','Labels']) #['Eke','Labels'])
Using default sph betas.
Summing over dims: {'it'}
Plotting from self.data[orb8][MFBLM], facetDims=['Eke', 'Labels'], pType=a with backend=pl.
Set plot to self.data['orb8']['plots']['MFBLM']['polar']
orbKey = 'orb7'
data.padPlot(keys=orbKey, Erange = [10, 15, 5], selDims={'Labels':tRot.round(2)[[0,2,5,7,9]]}, dataType='MFBLM', backend='pl', returnFlag=True,
height=400, width=1200, rc=[1,5]) #, facetDims=['Eke','Labels']) #['Eke','Labels'])
Using default sph betas. Summing over dims: set() Plotting from self.data[orb7][MFBLM], facetDims=['Eke', 'Labels'], pType=a with backend=pl. *** WARNING: plot dataset has min value < 0, min = (-0.0002881296903635873+3.586969769800971e-18j). This may be unphysical and/or result in plotting issues.
Set plot to self.data['orb7']['plots']['MFBLM']['polar']
orbKey = 'orb6'
data.padPlot(keys=orbKey, Erange = [10, 15, 5], selDims={'Labels':tRot.round(2)[[0,2,5,7,9]]}, dataType='MFBLM', backend='pl', returnFlag=True,
height=400, width=1200, rc=[1,5]) #, facetDims=['Eke','Labels']) #['Eke','Labels'])
Using default sph betas. Summing over dims: set() Plotting from self.data[orb6][MFBLM], facetDims=['Eke', 'Labels'], pType=a with backend=pl. *** WARNING: plot dataset has min value < 0, min = (-1.4300626655083803e-16-1.4065322861916882e-18j). This may be unphysical and/or result in plotting issues.
Set plot to self.data['orb6']['plots']['MFBLM']['polar']
These are from ePolyScat's getCro function, and are LF (unaligned ensemble) results. This provides a good, if general, overview.
# NEED TO SET AGAIN AFTER CLASS IMPORT!
import warnings
# warnings.filterwarnings('once') # Skip repeated numpy deprecation warnings in current build (xr15 env)
warnings.filterwarnings('ignore') # Skip repeated numpy deprecation warnings in current build (xr15 env)
# Comparative plot over datasets (all symmetries only)
Etype = 'Eke' # Set for Eke or Ehv energy scale
pGauge = 'L'
# Erange=[0, 100] # Plot range (full range if not passed to function below)
# data.plotGetCroComp(pType='SIGMA', Etype = Etype, Erange = Erange, backend = 'hv')
data.plotGetCroComp(pType='SIGMA', pGauge = pGauge, Etype = Etype)
# Comparative plot over datasets (all symmetries only)
data.plotGetCroComp(pType='BETA', Etype=Etype)
Same results as above, Ehv scale.
# Comparative plot over datasets (all symmetries only)
Etype = 'Ehv' # Set for Eke or Ehv energy scale
pGauge = 'L'
# Erange=[0, 100] # Plot range (full range if not passed to function below)
# data.plotGetCroComp(pType='SIGMA', Etype = Etype, Erange = Erange, backend = 'hv')
data.plotGetCroComp(pType='SIGMA', pGauge = pGauge, Etype = Etype)
# Comparative plot over datasets (all symmetries only)
data.plotGetCroComp(pType='BETA', Etype=Etype)
# TODO
# data.plotGetCroComp(pType='BETA', Etype=Etype, backend='hv')
# Betas vs. Gauge
data.plotGetCro(pType='BETA', Etype=Etype)
# print(data.jobInfo['ePolyScat'][0])
# print('Run: ' + jobInfo['Starting'][0].split('at')[1])
import scooby
scooby.Report(additional=['epsproc', 'xarray', 'jupyter'])
| Thu Aug 29 12:14:21 2024 EDT | |||||||
| OS | Linux | CPU(s) | 64 | Machine | x86_64 | Architecture | 64bit |
| RAM | 62.8 GiB | Environment | Jupyter | File system | btrfs | ||
| Python 3.10.11 | packaged by conda-forge | (main, May 10 2023, 18:58:44) [GCC 11.3.0] | |||||||
| epsproc | 1.3.2-dev | xarray | 2022.3.0 | jupyter | Version unknown | numpy | 1.23.5 |
| scipy | 1.10.1 | IPython | 8.13.2 | matplotlib | 3.5.3 | scooby | 0.7.2 |
# Check current Git commit for local ePSproc version
!git -C {Path(ep.__file__).parent} branch
!git -C {Path(ep.__file__).parent} log --format="%H" -n 1
* 3d-AFPAD-dev
dev
master
93fc0fb4520e1bd89cfebba8f5354b6791c37a1d
# Check current remote commits
!git ls-remote --heads https://github.com/phockett/ePSproc
8c9201fdcb98f42875730b94282f50328decab8a refs/heads/3d-AFPAD-dev 7e4270370d66df44c334675ac487c87d702408da refs/heads/dev 1c0b8fd409648f07c85f4f20628b5ea7627e0c4e refs/heads/master 69cd89ce5bc0ad6d465a4bd8df6fba15d3fd1aee refs/heads/numba-tests ea30878c842f09d525fbf39fa269fa2302a13b57 refs/heads/revert-9-master baf0be0c962e8ab3c3df57c8f70f0e939f99cbd7 refs/heads/testDev
!hostname
8fdb9101787f
!conda env list
# conda environments: # base * /opt/conda